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Computational  efficiency  is  highly  important  for  upscaling  detailed  electrode-level  and  cell-level  models 
to  the  system  level  required  for  the  design  and  control  of  fuel  cells.  We  present  a  computationally  efficient 
1 D  + 1 D  fuel  cell  model  based  on  a  combination  of  analytical  and  numerical  approaches.  On  the  electrode 
level,  we  develop  approximate  analytical  solutions  for  the  ID  current/potential  distribution  via  a  hybrid 
algorithm  of  power-law  approach  and  perturbation  method.  Compared  to  the  conventional  perturbation 
method,  this  work  keeps  the  intrinsic  nonlinearity  of  electrochemical  kinetics,  while  providing  clearer 
physical  meaning  than  some  purely  mathematical  methods  like  the  Adomian  decomposition  method.  By 
integrating  the  resulting  overpotential  profile  into  mass  transfer  models,  concentration  overpotentials 
are  obtained  and  the  thermodynamic  framework  is  then  used  for  analyzing  the  H2/CO  electrochemical 
co-oxidation  kinetics.  A  novel  expression  is  also  presented  to  interconvert  volume-  and  area-specific 
exchange  current  densities.  On  the  cell  level,  a  linear  relationship  between  local  current  density  and  solid 
temperature  is  further  developed  for  efficient  ID  +  ID  thermal  along-the-channel  numerical  simulations 
without  requiring  computational  iterations.  Both  the  electrode-level  and  cell-level  macroscopic  fuel  cell 
models  are  validated  against  full  numerical  solutions  available  in  the  open  literatures  over  a  wide  range 
of  operating  conditions.  With  the  hybrid  analytical/numerical  approximation  in  two  dimensions,  the 
computational  framework  is  predicted  to  be  sufficiently  efficient  for  real-time  simulations. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  environmentally  benign,  high-efficiency  energy 
conversion  devices,  and  have  been  considered  to  be  potential  candi¬ 
dates  for  stationary  power  plants,  automotive  engines  and  portable 
applications.  Modeling  and  simulation  techniques  are  becoming 
increasingly  important  for  both,  understanding  fundamental  pro¬ 
cesses  in  fuel  cells,  and  supporting  design  and  optimization  of  fuel 
cell  systems. 

In  recent  years,  many  researchers  have  been  developing 
multi-scale  modeling  frameworks  of  fuel  cells  [1-5].  Microscopic 
modeling  aims  at  a  physics-oriented  explanation  of  the  ele¬ 
mentary  kinetics  of  electrochemistry  and  transport,  including  its 
impedance-based  analysis  [6-8].  Macroscopic  models  have  been 
developed  for  multi-dimensional,  non-isothermal  and  transient 
simulation  with  computational  fluid  dynamics  (CFD)  based  codes 
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[9-11].  Although  a  direct  combination  of  microscopic  and  macro¬ 
scopic  models  would  provide  a  better  fundamental  understanding 
of  fuel  cells,  the  complexity  of  both  approaches  makes  it  generally 
difficult  to  couple  them  into  one  single  multi-scale  model  due  to 
high  numerical  cost.  On  the  other  hand,  simplified  semi-empirical 
and  analytical  models  are  widely  used  in  predicting  polarization 
behavior  of  fuel  cells  [12,13].  Flowever,  the  coarse  framework  and 
lacking  of  physical  meaning  limits  their  application,  especially  for 
spatially  distributed  modeling. 

An  accurate  electrode-level  and  cell-level  model  with  low 
computational  cost  is  essential  for  system-level  analysis,  espe¬ 
cially  for  controller  design  and  hardware-in-the-loop  simulations. 
Approximate  analytical  solutions  provide  a  good  balance  between 
mechanistic  and  semi-empirical  models.  Gurau  [14]  presented  an 
analytical  solution  of  a  transport  model  for  polymer  electrolyte 
membrane  fuel  cells  (PEMFCs)  with  the  assumption  of  constant 
overpotential  in  the  catalyst  layer.  Considering  the  catalyst  layer  as 
an  infinitely  thin  interface,  Tsai  [15]  proposed  a  two-dimensional 
analytical  expression  of  oxygen  transport  in  the  cathode  diffusion 
layer  of  a  PEMFC.  Costamagna  [16]  obtained  an  analytical  solu¬ 
tion  of  overpotential  distribution  in  solid  oxide  fuel  cell  (SOFC) 
electrodes  by  assuming  constant  gas  concentration  and  exchange 
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Nomenclature 

A  area(m2) 

[A],[B]  matrixes  of  linear  relationship,  [I]  =  [/V ] [ Ts  ]  +  [B] 

C  concentration  (mol  nrr3 ) 

Cp  molar  specific  heat  (J  mol-1  K-1 ) 

D  diffusivity  (m2  s)  or  channel  depth  (m) 

F  Faraday’s  constant  (96,487  C  mol-1 ) 

H  molar  enthalpy  (J  mol-1 ) 

z'o  volume-specific  exchange  current  density  (A  nrr3 ) 

/0  area-specific  exchange  current  density  (A  nrr 2 ) 

/  current  density  (A  nrr 2 ) 

j  electrochemical  reaction  rate  (A  nr 3 ) 

/<o  dimensionless  variable  in  Eq.  (3) 

l  electrode  thickness  (m) 

L  cell  length  (m) 

n  power-law  index  or  number  of  control  volumes 

rze  electrons  transferred  per  reacting  molecule,  ne  =  2 

N  flux  (mol  nr2  s-1) 

p  pressure  (Pa) 

r  volumetric  reaction  rate  (mol  nr3  s-1 ) 

R  universal  gas  constant  (8.314Jmol_1  K-1)  or  reac¬ 

tion  rate  per  area  (mol  nrr2  s-1 ) 

S  molar  entropy  (J  mol-1  K-1 ) 

T  temperature  (K) 

u  velocity  (m  s-1 ) 

W  width  (m) 

x  molar  fraction  or  coordinate  in  thickness  direction 

(m) 

y  dimensionless  overpotential,  y  =  neFr]/RT 

z  coordinate  in  the  cell  length  thickness  (m) 

Greek  letters 

aa  anodic  transfer  coefficient 

ac  cathodic  transfer  coefficient 

r)  overpotential  (V) 

k  heat  conductivity  (W  m-1  K-1 ) 

0  potential  (V) 

p  density  (kg nr3) 

a  conductivity  (S  nr1 ) 

v  stoichiometric  coefficient  of  reaction 

Subscripts  and  superscripts 
a  anode 

c  cathode 

e  electrolyte 

eff  effective 

el  electronic  conducting  phase 

h  heat  transfer 

i,  j  species 

in  inlet 

ion  ionic  conducting  phase 

m  mass  transfer 

ref  reference  or  reforming  reaction 

s  solid  phase 

WGS  water  gas  shift  reaction 


current  density.  Considering  variation  of  both  gas  concentration 
and  overpotential,  Bao  [17]  further  proposed  an  approximate  ana¬ 
lytical  solution  of  a  transport  model  for  anode-supported  SOFCs. 
Other  approximate  analytical  expressions  were  also  developed  for 
modeling  channel  flow  and  heat  transfer  in  a  fuel  cell  [18,19].  An 
overview  of  analytical  fuel  cell  models  is  given  in  Kulikovsky’s 
textbook  [20]. 


This  paper  will  firstly  introduce  a  continuum  transport  model 
in  porous  electrodes  of  a  fuel  cell.  Depending  on  the  electrode 
thickness,  different  approximate  analytical  solutions  of  overpoten¬ 
tial  distribution,  charge  transfer,  and  electrochemical  reaction  are 
developed  and  then  applied  for  modeling  the  mass  transfer  for 
various  reactant  systems.  Finally,  the  ID  electrode  model  will  be 
introduced  into  a  numerical  1 D  along- the-channel  cell-level  model. 
An  efficient  computational  framework  is  developed  between  the 
local  current  density  and  cell  temperature.  The  models  and  results 
are  discussed  in  the  context  of  SOFCs,  and  are  partly  suitable  for 
PEMFCs. 


2.  Background:  electrode-level  charge  transport  model 


The  core  structure  of  a  fuel  cell  is  a  “sandwich”  configuration 
of  the  two  electrodes  and  the  membrane  (membrane-electrode 
assembly,  MEA).  The  electrodes  of  an  SOFC,  or  the  catalyst  layers 
of  a  PEMFC,  are  formed  by  a  mixture  of  an  ionic  (ion)  conductor, 
an  electronic  (el)  conductor,  and  pore  space  for  gas  diffusion,  the 
connection  of  which  make  up  three-phase  boundary  (TPB)  regions. 
Close  to  the  TPB,  electrochemical  reactions  occur,  exchanging 
electrical  charges  between  ionic  conducting  phase  and  electronic 
conducting  phase. 

The  model  developed  here  is  based  on  the  following  assump¬ 
tions:  the  gas-phase  is  ideal,  temperature  and  pressure  are  uniform 
throughout  the  electrode,  both  of  the  two  conducting  phases  are 
continuous  and  homogeneous,  and  only  steady  state  is  considered. 
Furthermore,  in  order  to  decouple  the  problems  of  charge  trans¬ 
fer  and  mass  transfer  for  computational  efficiency,  we  neglect  the 
influence  of  gas  concentration  on  the  exchange  current  density  and 
electrochemical  reaction. 

Along  the  positive  direction  from  electrode/channel  (E/C) 
interface  to  electrode/electrolyte  (E/E)  interface  (x  coordinate), 
according  to  Ohm’s  law,  the  charge  transfer  in  two  phases  is 


^  '  (  ^ion^^ion)  —  ^  '  (^el^^el)  —  j 


(1) 


where  the  0  and  a  are  the  potential  and  conductivity,  respec¬ 
tively,  and  the  electrochemical  reaction  rate  j  can  be  described  by 
a  Butler-Volmer  (BV)  equation, 


J  =  i  o 


(aaneF  \ 

exp(-Try- 


'0  =  'O.ref  exp 


£act 

~R~ 


exp  [  - 
1 

T  4  f 


(  ac neF  \ 
(  RT  7 


RT 


(2) 


where  F  is  the  Faraday  constant,  R  is  the  universal  gas  constant,  T 
is  the  operating  temperature,  ne  is  the  number  of  electrons  partic¬ 
ipating  in  the  electrochemical  reaction,  aa  and  ac  are  the  anodic 
and  cathodic  charge  transfer  coefficients,  z'o  is  the  volume-specific 
exchange  current  density  (which  depends  on  both  the  TPB-specific 
exchange  current  density  and  the  microstructure  of  the  electrode), 
z'q  ref  is  the  reference  exchange  current  density  at  the  reference 
temperature  Tref,  Eact  is  the  activation  energy,  p*  and  v,-  are  the 
partial  pressure  and  reaction  order  of  species  z.  The  overpotential 
r\  is  defined  as  the  potential  difference  between  the  two  phases 
(0ei  -  0jon)  minus  that  in  equilibrium.  At  the  E/C  interface,  the  ionic 
current  density  completely  transfers  into  the  electronic  current 
density,  while  at  the  E/E  interface,  the  electronic  current  density 
completely  transfers  into  the  ionic  current  density. 

The  exchange  current  density  z0  depends  on  species  concentra¬ 
tion  (last  term  in  Eq.  (2)),  which  is  generally  a  function  of  position 
in  the  electrode  thickness  (here:  x  dimension)  and  along  channel 
length  (here:  z  dimension).  In  the  present  work,  in  order  to  derive 
computationally  efficient  analytical  solutions,  we  assume  that  z0 
is  constant  through  the  electrode  thickness  and  varies  only  along 
the  channel  length:  the  species  bulk  partial  pressure  at  E/C  interface 
plib(z)  rather  than  the  TPB  concentrations  p,(x,z)  is  used  in  Eq.  (2)  for 
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the  concentration-dependent  exchange  current  density.  Note  that 
this  assumption  only  influences  the  activation  overpotential.  The 
concentration  overpotential  is  calculated  using  the  fully  resolved 
concentrations  (Eq.  (59)). 

Define  the  following  dimensionless  variables 

f=y,  y=fn .  *b  =  ^  0) 

l  crt 

where /=  neF/RT,  crt  =  crion<jeil(<7ion  +  crel),  l  is  the  electrode  thickness, 
£  and  y  are  dimensionless  electrode  thickness  and  overpotential, 
respectively,  we  obtain  from  Eqs.  (1 )  and  (2)  the  governing  equation 
for  overpotential  distribution  in  a  composite  electrode, 


y  =  k0[exp(aay)  -  exp(-acy)] 

/(?  =  0)  =  N0  =  -jr-,  y(?=l)  =  Ni  =  ^- 


(4) 


where  I  is  the  operating  current  density. 

There  are  two  main  difficulties  to  solve  this  system  analytically, 
one  lies  in  the  intrinsic  equation  structure  related  to  the  variable 
/<o,  the  other  is  the  strong  nonlinearity  of  exponent  function.  In 
addition,  lacking  of  Dirichlet-type  boundary  conditions  also  brings 
some  difficulties  to  this  boundary-value  problem,  which  will  be 
explained  later. 

The  variable  k0  is  similar  to  the  Thiele  modulus  for  a  classical 
problem  of  diffusion  and  reaction  of  gases  throughout  a  porous  cat¬ 
alyst  pellet  in  isothermal  conditions.  In  the  present  situation,  low 
k0  values  mean  that  the  electrode  performance  is  kinetics-limited, 
whereas  high  k0  values  mean  that  the  electrode  performance  is  lim¬ 
ited  by  the  ionic  conductivity.  When  k0  is  small  enough  (k0  ->  0), 
the  system  reduces  to  a  linear  one,  which  does  not  satisfy  the  two 
Newman-type  boundaries  at  the  same  time.  When  k0  is  big  enough 
(1/fco  -*  0),  by  multiplying  l/k0  at  both  sides  of  Eq.  (4),  the  differ¬ 
ential  system  reduces  to  an  algebraic  one,  which  means  a  strong 
boundary  effect;  this  has  been  discussed  in  our  previous  work  [17]. 
Because  k0  is  proportional  to  l 2,  both  cases  occur  in  practice  due  to 
the  wide  range  of  electrode  size,  e.g.  the  anode  in  anode-supported 
SOFC  and  the  catalyst  layer  in  PEMFC. 

In  literature  [16,17,20],  an  exact  solution  of  Eq.  (4)  has  typically 
been  obtained  via  linearization  of  the  exponent,  that  is,  by  assuming 
that  overpotentials  are  low.  It  results  in 


yiin  — 


Ni-Noe-kM  .  Nj-Npe*-  u 
X(ex-e~x)  X(ex-e~x) 


X  —  \J ko(oia  +  oic )  (5) 


Based  on  the  singular  perturbation  method  (SPM),  we  have  pre¬ 
viously  obtained  an  approximate  profile  of  overpotential,  which 
consists  of  a  logarithmic  term  due  to  mass  transfer  and  two  expo¬ 
nent  terms  due  to  boundary  effects  in  a  thick  electrode  [17].  When 
the  variation  of  gas  concentration  is  neglected  (the  logarithmic 
term  vanishes),  it  reduces  to 


N  i 

yspm  —  -j-e 


-Mi-0 


(6) 


Although  this  expression  provides  a  relatively  compact  form 
and  clearer  physical  meaning,  it  is  almost  numerically  identical  to 
Eq.  (5)  when  e~x  ->  0.  Due  to  the  inherent  approximation  of  lin¬ 
earization  the  exponent  term,  both  of  them  are  inaccurate  when 
the  overpotential  is  not  small.  The  key  problem  in  this  case  is  how 
to  keep  the  nonlinearity  of  system. 


operators  [21  ].  For  convenience  to  the  reader,  we  briefly  introduce 
the  principle  of  ADM.  Consider  the  generalized  differential  equation 

Ly  +  Ry  +  Ny  =  g(0  (7) 

where  L  is  the  highest  order  derivative,  Ry  is  the  linear  differential 
operator  of  less  order  than  L,  Ny  represents  the  nonlinear  terms, 
and  g(£ )  is  the  source  term.  Applying  the  inverse  operator  L-1  to 
both  sides  of  Eq.  (7),  we  obtain 

y=m)-L-\Ry)-L-\Ny )  (8) 

where  the  function/(£)  represents  the  terms  arising  from  integrat¬ 
ing  g(f).  Next,  the  nonlinear  operator  Ny  =  F(y)  is  represented  by  an 
infinite  series  of  the  so-called  Adomian  polynomials 


f(y)  =  YjA”  O) 

n= 0 

Define  the  solution  y(f)  by  the  series 


oo 

y  =  y /n  (io) 

n= 0 

where  the  components  yn  are  determined  recursively  by  the  mod¬ 
ified  ADM  [22] 


yo  =/o(f)> 

yi  =MS)-L-\Ry0)-L-\Ao),  (ID 

yu+2  =  -L~HRyk+\)-L~\Ak+-i),  k> o 


where  /(£)=/o(0+/i(0>  and  the  Adomian  polynomials  are  con¬ 
structed  by  the  following  algorithm 


Myo,y\, 


(12) 


In  the  fuel  cell  electrode  model  presented  above  (Eq.  (4)),  L  is  the 
2nd-order  differential  and  L-1  the  two-fold  integral,  Ry  =  0,  g(ij)  =  0 
and  Ny  =  -k0[ exp ( aay)  -  exp(-o'cy)].  Therefore, 


yo  =  c0, 

yi  =  -— f  +  ^(e“aC°  -  e-“cC°)£2, 
Oel  2 


y2  =  ^(aaeaaC°  +ace  acC°) 
6 


-MLf  +  ^(e^co  _  e-«cc0)^4 
L  ae\  4 


(13) 


Substituting  into  Eq.  ( 1 0),  we  obtain  the  ADM  solution  of  overpo¬ 
tential,  where  the  constant  Co  can  be  determined  by  the  boundary 
condition,  dy/df|^_1  =N^.  Fig.  1  shows  the  comparison  between 
the  exact  solution  and  ADM  approximation  with  a  different  number 
of  terms.  With  more  terms,  the  ADM  solution  approaches  gradually 
to  the  exact  solution.  When  k0  is  not  large,  e.g.  for  the  catalyst  layer 
in  PEMFC,  ADM  is  an  effective  method  for  approximation.  How¬ 
ever,  more  terms  are  generally  required  to  ensure  accuracy,  which 
leads  a  complicated  expression  without  clear  physical  meaning. 
Furthermore,  it  is  hard  to  obtain  an  explicit  expression  due  to  the 
implicit  iteration  of  constant  c0.  Under  the  condition  of  a  large  k0 
(/< o  >200),  for  example,  for  anode-supported  SOFCs,  we  found  the 
ADM  solution  not  satisfying  and  even  not  convergent. 


3.  Nonlinear  approximation  of  the  electrode-level  model 

3.1.  Adomian  decomposition  method 

The  Adomian  decomposition  method  (ADM)  allows  solution  of 
nonlinear  functional  equations  of  various  kinds  (algebraic,  differ¬ 
ential,  partial  differential,  integral,  etc.)  without  approximating  the 


3.2.  Power-law  approach 

For  thick  electrodes,  as  present  for  example  in  anode-supported 
SOFCs,  1  /k0  is  small  enough  to  be  a  perturbation  variable  [  1 7 ],  which 
can  be  applied  to  simplify  the  problem.  Here,  we  employ  a  power 
function  to  approximate  the  exponential  functions,  that  is, 

e«a y  -e-^y  =  k{y  +  a)n 


(14) 
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Dimensionless  electrode  thickness 


Fig.  1.  Comparison  between  the  exact  solution  and  ADM  solution  with  a  different 
number  of  terms  for  parameters  corresponding  to  an  intermediate-temperature 
SOFC:  /  =  50 |jim,  aa  =  ac  =  0.5,  i0  =  5  x  10s Arrr3,  1= 4xl04Arrr2,  7=1073.15I<, 
cr;on  =  2.267  S  m_\  <xei  =  3.03  x  1 04  S  rrr 1  and  ko  =  1 1 .93. 


where  k,  a,  n  are  parameters  to  be  determined  later.  Substituting 
into  Eq.  (4)  yields 

/  y"  =  k0k{y  +  a)n  t=y+a  f  t"  =  k0ktn  (  , 

\y(l)  =  N!  \t,(l)  =  N1  1  j 


The  solution  of  the  above  system  is  a  power-law  expression, 
that  is,  t  =  t(l)[l  +b{  1  -£)]m.  Substituting  into  Eq.  (15),  we  obtain 
m  =  2/(l  -  n)  and 


k 


/<0(1  -  ny 


[y(i)  +  a]n 


2(1  +n) 

-2b[y(l)  +  a]  =  (1  -  n)Ni 


(16) 


In  order  to  obtain  values  for  the  coefficients,  both  the  zero-order 
and  first-order  differential  of  the  power  function  approach  is  set  to 
that  of  the  corresponding  full  exponential  expression  at  the  E/E 
interface. 

h  =e“^(D-e-“^(1)  =  k[y(l)  +  a]n 

h  =  aiae^O)  +  ace-a =  nk[y(  1)  +  a]"-1  1  j 


The  explicit  expressions  of  n,  a,  k,  b  can  be  solved  from  Eqs.  (16) 
and  (17) 


k  = 


hN[ 

-/2JV2  ’ 

/l 


/iN? 


2ko/,2  -/2N? 


Cf1N12/(2kQf12-/2N12))n’ 


b  =  - 


2N1(2k0/12  -/2N2) 
(l-nF,N2 


(18) 


In  order  to  determine  y(l),  we  reenter  into  the  original  system 
of  Eq.  (4)  to  obtain  an  exact  implicit  solution 


(y')2  =  2k0(l< 


1 


'  +  C 


(19) 


As  mentioned  above,  lacking  Dirichlet-type  boundary  condition 
brings  difficulty  to  determine  the  constant  C.  Considering  the  sign  of 
1  st-  and  2nd-order  derivatives  at  the  electrode  boundaries  ( d2y  >  0, 
dy(0)  <  0,  dy(l )  >  0),  the  solution  should  show  a  concave  profile.  For 
a  thick  electrode,  most  of  the  electrode  is  inactive  to  electrochem¬ 
ical  reaction,  which  means  at  the  minimum  point 


(20) 


Therefore,  C=  -(l/aa  +  l/ac).  and  the  overpotential  at  the  E/E 
interface  y(l )  can  be  solved  by 


2 k0  ( 3-e“ay(i)  +  JLe-«cy(U  _  JL  _  JA  =  n2 

V  Ola.  Oic  Ota  OtcJ  1 

When  eta  =  etc,  there  is  an  exact  solution 


1 


y(D  =  ->nU 


2k0 


+  2  + 


(21) 


(22) 


Thus,  we  get  an  explicit  power  function  which  is  valid  in  the 
region  close  to  the  E/E  interface 

yPi  =  [y(l)  +  a][l  +b(l  -4T)]2/(1-n)  -a  (23) 

The  power  function  close  to  the  E/C  boundary  can  be  obtained 
similarly  and  is  given  by 

yPo  =  [y(0)  +  a'][l  -  b'£]2/(1“n,)  -  a'  (24) 

where  a',  b',  n',y( 0)  are  calculated  by  replacing  Nlty(l)  to  N0,y(0) 
in  all  the  related  equations. 


3.3.  Hybrid  power  law  and  linear  approximation 


The  power  law  approximation  as  derived  above  is  valid  only  in 
the  region  close  the  electrode  interfaces.  On  the  other  hand,  the  lin¬ 
ear  approximation  (Eq.  (5))  should  be  valid  in  most  of  the  region  in  a 
thick  electrode  due  to  a  small  electrochemical  reaction  rate.  There¬ 
fore,  we  propose  a  hybrid  function  to  make  a  smooth  transition 
between  these  two  expressions. 

Based  on  the  analysis  of  boundary-layer  effect  in  our  previous 
work  [17],  we  construct 

y  =yiin,WKB  +(yPl  -yiin,WKB)e_S°A(1_^  +&l(0  (25) 


where 


yiin,WKB 


(Nt  -  N0e~sox)eso^  +  (Ni  -  N0es o^)e~so^ 
SoA,(eso^  -  e_so^) 


s0  = 


eaa  _  e-ac 
Ota  "T  Ot c 


(26) 


Here,  we  apply  an  additional  correction  to  the  linear  approxima¬ 
tion  (Eq.  (5))  based  on  the  WKB  (Wentzel,  Kramers  and  Brillouin) 
perturbation  method  [23],  which  is  explained  in  detail  in  Appendix 
A.  In  order  to  make  y  strictly  satisfy  the  boundary  conditions  in  Eq. 
(4),  gi(0  is  added  and  derived  by  considering  the  requirements 

gi(  0)  =  0,  si(0)  =  0,  gi(l)  =  0, 

g'lW  =  -So^tVpl(l)  -yiin,WKB(^)l  (27) 

yielding 

gi(4T)  =  s0X[ypl(l) -ylin>WI<B(l)]4T(l  _£)e-siMW)  (28) 

In  principle,  the  slope  rate  Si  should  be  fitted  from  the  exact 
numerical  solution.  However,  we  tried  to  simply  chose  Si  =  1.3so, 
which  is  satisfying  enough  in  a  wide  range  of  applications. 

The  final  approximation  solution  using  the  hybrid  power-law 
and  linearization  approach  is 

y  =yiin,WKB  +(yPl  -yiin,WKB)e_S°^1_^  +  S(A[yPl  (1  )  -  yiin,WKB(l )] 
?(1  _  f)e-Sl^1-^  +  (yp0  —  yiin,WKB)^_s°^^  +  so-MyPo(0) 
-yiin,WKB(°)]^(^  -  Oe~S^  (29) 


In  practical  fuel  cells,  the  electronic  conductivity  is  generally 
much  larger  than  the  ionic  conductivity  (crion  <£0^1  )*  the  boundary 


y'  =  0, 


y^0 
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Dimensionless  electrode  thickness  Dimensionless  electrode  thickness 


Fig.  2.  Comparison  between  the  exact  solution  and  different  approximation  solutions  for  parameters  corresponding  to  an  anode-supported  SOFC:  1=1  mm,  aa  =  ac  =  0.5, 
io  =  5  x  108  Am-3, 1  =  4  x  104Arrr2,  7=1073. 15K,  crion  =  2.267  Sirr1,  crel=3.03  x  104Sm-1  and  /<o=4771.  (a)  Linear  presentation  of  overpotential  as  function  of  electrode 
thickness,  (b)  logarithmic  presentation  of  the  same  data,  and  (c)  linear  presentation  of  the  region  close  to  E/E  interface. 


effect  at  the  E/C  interface  is  negligible  and  only  the  first  three  terms 
in  Eq.  (29)  are  needed. 

For  a  typical  anode-supported  SOFC,  Fig.  2(a)  shows  the  overpo¬ 
tential  profiles  calculated  via  the  exact  numerical  solution,  linear 
approximation  in  Eq.  (5),  SPM  approximation  in  Eq.  (6),  WKB 
approximation  in  Eq.  (26),  hybrid  power-law  approximation  in 


Fig.  3.  Index  of  the  power  function  (Eq.  (14))  as  function  of  the  operating  current 
density,  other  parameters  are  identical  to  those  in  Fig.  2. 


Eq.  (29)  and  the  power  function  in  Eq.  (23).  All  the  overpoten¬ 
tial  profiles  show  a  strong  boundary-layer  effect  due  to  a  large  k0. 
For  a  clearer  image,  Fig.  2(b)  and  (c)  further  shows  the  logarith¬ 
mic  or  linear  overpotential  profile  within  the  whole  electrode  or 
the  region  close  to  the  E/E  interface.  The  linear  and  SPM  approx¬ 
imations  are  almost  overlapping  as  mentioned  before.  Note  that 
although  starting  from  the  assumption  of  a  thick  electrode,  the 
power-law  approach  is  found  to  be  accurate  as  k0  >  20,  which  cor¬ 
responds  to  an  electrode  thickness  of  around  60  pun  assuming 
other  typical  parameters.  Fig.  3  shows  the  effect  of  the  operat¬ 
ing  current  density  on  the  index  of  power  function.  As  current 
density  decreases,  the  power-law  index  approaches  gradually  to 
unity,  which  means  a  linear  approximation.  Fig.  4(a)  and  (b)  respec¬ 
tively  shows  the  linear-  and  logarithmic-type  exact  and  power-law 
solutions  when  the  electronic  conductivity  is  equal  to  the  ionic  con¬ 
ductivity  (crion  =  crej).  Such  a  situation  is  present  in  the  so-called 
IDEAL-cell  concept  [24],  as  a  result,  all  the  five  terms  in  Eq.  (29)  are 
necessary,  thus  the  boundary-layer  effects  at  both  sides  of  electrode 
can  be  accurately  reflected  by  the  power-law  approximation.  Due 
to  analytical  expression,  the  hybrid  power-law  and  linear  approx¬ 
imation  is  smoother  than  the  exact  numeric  solution  with  250 
finite-difference  discrete  grids. 

3.4.  Relationship  between  volume-specific  and  area-specific 
exchange  current  densities 

In  the  above  electrode  model,  electrochemical  reactions  are 
assumed  to  occur  throughout  the  complete  thickness  of  the  porous 
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Fig.  4.  Comparison  between  the  exact  solution  and  power-law  approximation  for 
crei  =  cr and  1  =  3  x  104  Am-2,  other  parameters  are  identical  to  those  in  Fig.  2.  (a) 
Linear  and  (b)  logarithmic  presentation. 


electrode  (volumetric  model).  The  exchange  current  density  io  is 
a  volume-specific  parameter  (given  in  Ampere  per  electrode  vol¬ 
ume).  On  the  other  hand,  many  models  in  literature  assume  that 
electrochemistry  takes  only  place  at  the  E/E  interface  (interfacial 
model).  The  exchange  current  density  used  in  those  models  is  an 
area-specific  parameter  (given  in  Ampere  per  cell  area),  that  is, 

I  =  JoIexpCo-aftact)  -  exp(-a</i?act)]  (30) 


where  rjact  is  the  total  activation  overpotential,  /0  the  area-specific 
exchange  current  density. 

Eq.  (19)  provides  an  exact  solution  of  the  overpotential  at  the 
E/E  interface,  r)EfE.  As  crion<<rel,  ?7e/e  is  almost  equal  to  the  total 
activation  overpotential  (refer  to  Eq.  (59)).  Due  to  equivalence  of 
these  two  modeling  approaches,  r)act  =  77e/e>  a  relationship  between 
volume-specific  and  area-specific  exchange  current  densities  can 
be  derived.  For  aa  =  ac , 


_ Qjfft/f2 _ 

2<(\/W  +  4  +  aac) 


(31) 


For  a  thick  electrode  (Eq.  (20)  is  valid),  aaC=-2.  Expand  the 
square  root  term  based  on  the  binomial  theorem  as  generally 


(J/J0)2  « 1,  we  obtain  an  expression  independent  of  the  operation 
current  density 

h_2J^olx2Jlal  (32) 

^ion  aion 

This  expression  shows  a  nonlinear  relationship  between  these 
two  exchange  current  densities,  which  is  different  to  the  common 
relationship  with  a  simple  geometric  factor,  i’o  =  /o/geo,  which  we 
believe  lacks  theoretical  support. 


4.  Integration  of  mass  transfer  into  electrode-level  model 


The  model  presented  in  the  sections  above  treats  charge  transfer 
only.  Here,  we  add  expressions  for  mass  transfer.  Note  that,  in  order 
to  decouple  mass  and  charge  transport  formulations,  the  influence 
of  concentration  on  exchange  current  density  is  neglected  in  this 
work.  The  mass  transfer  of  multi-component  gas  species  in  a  porous 
electrode  is  governed  by  the  Stefan-Maxwell  equations, 


V*.  =  E 


xtNj  -  XjNt 
cDjjeff 


(33) 


where  c  is  the  concentration  of  gas  mixture,  x*  =  q/ct  and  Nf  are 
the  molar  fraction  and  flux  of  species  i,  respectively,  Dy  eff  is  the 
effective  diffusivity,  which  includes  the  binary  diffusivity  between 
species  i  and  j  (Dy)  and  their  Knudsen  diffusion  coefficients  (DlK, 
DjK)  and  is  corrected  by  the  electrode  porosity  (e)  and  tortuosity 
Was  [2] 


Oy.eff  — 


2r2 


(1/D,j)  +  (1/D,-k)  (1/Dy)  +  (1/Dj  K) 


(34) 


where  rP  is  the  average  particle  size  of  electrode  and  M*  is  the 
molecular  weight  of  species  i. 

The  mass  balance  of  species  i  is  given  by 


v-(n«)=0+E^  (35) 

,lc  k 

where  iq,  uI  k  are  the  stoichiometric  coefficients  of  species  i  in  the 
electrochemical  reaction  and  in  chemical  reaction  k  (e.g.,  reforming 
reaction  and  water  gas  shift  reaction),  respectively. 

As  boundary  conditions,  the  species  concentration  at  the  E/C 
interface  is  assumed  to  be  equal  to  the  bulk  channel  concentration, 
and  at  the  E/E  interface,  the  species  can  not  penetrate  the  dense 
electrolyte  layer  (note  this  assumption  is  not  correct  for  PEMFCs, 
where  there  is  complex  water  management  in  the  PEM),  that  is, 

xi(x  =  0)  =  Xjb,  Ni(x  =  /)  =  0. 


4.1.  H2-H20  or  C0-C02  system 


In  this  binary  stationary  system  (1:  H2,  2:  H20),  equimolar 
counter  diffusion  occurs  ( N-i=-N2 )  and  there  are  no  chemical 
reactions,  that  is,  dN^ldx  =  v^jlneF  or  d2Xi/dx2  = -Uij7cD12effneF. 
Because  d2r//dx2  =jlcrt}  there  is  a  simple  relationship  between  the 
profiles  of  gas  concentration  and  overpotential 


*i(*)  =  *i,b- 


V\CTt 


c^12,effneF 


T]{X)-T]{0)-  — X 
Oinn 


(36) 


where  ;/(x)  has  been  determined  from  the  power-law  approach. 
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4.2.  H2-H2O-N2  or  C0-C02-N2  system 

In  this  case  (1:  H2,  2:  H20,  3:  N2),  N2  is  inert  (N3  =  0, 
note  that  we  investigate  the  stationary  case  only)  and  H2- 
H20  shows  an  equimolar  counter  diffusion  (N-i=-N2),  i.e. 
d2(lnx3)/dx2  =  L>ij(l/D13ieff-  l/D23,eff)/cncF.  So, 


where  D*  is  the  Fick  diffusivity  of  species  i,  which  is  considered 
constant  for  simplicity.  Because  Fields  law  is  not  strictly  valid 
for  multi-component  transport,  N2  is  treated  here  as  a  balancing 
species,  that  is,  the  above  equation  is  not  used  for  N2  and  the  N2 
molar  fraction  (x5 )  follows  from  Y^xi =  1  • 

Consider  electrochemical  oxidations  of  both  FI2  and  CO, 


x3(x) 


=  x3,b  exp 


f  VtV\ 
\cneF 


1 

^13,eff 


1 

^23,eff 


rj{x)  -  rj(0) - x 

*7jnn 


Xi(x) 

=  *l,b 


^23,eff[^13,eff  ln(x3A3,b)  +  (^12,eff  ~  ^13,eff)(x3  ~  *3,b)l 
^12,eff(^23,eff  “  ^13,eff) 


(38) 


When  x3b  =  0,  Eq.  (38)  reduces  to  Eq.  (36)  automatically. 
4.3.  02-N2  system 


dNj_  _  _dN2  _  _ajj_  dN3  _  dN4  _  (1  -  co)j 

dx  ~  dx  ~  neF+rwGS’  dx  ~  dx  ~  neF 

-  rwes  (44) 

where  00  is  the  proportion  of  FI2  electrochemical  current  density  in 
the  total  current  density,  00  =  Jh2/C/h2  +  jco).  and  rWGS  is  the  rate  of 
the  water-gas  shift  reaction  (C0  +  H20^C02  +H2),  which  is  con¬ 
sidered  as  the  only  chemical  reaction  in  this  paper.  Note  that  in  the 
present  work  we  neglect  carbon  formation  (e.g.,  due  to  Boudouard 
reaction);  also,  because  no  CH4  is  assumed  present,  there  are  no 
steam  reforming  or  dry  reforming  reactions. 

At  the  E/E  interface,  fluxes  of  all  the  species  are  zero 
(Nj|  =  0),  which  leads  to  Ni=-N2,  N3  =  -N4,  £Nf  =  0, 
dN !  /dx  +  dN3/dx  =  - cD !  d2Xi  /dx2  -  cD3d2x3 /dx2  =  -jlneF.  Therefore, 

DiXi  +  03X3  =  £f*-  s5bx + a =/(x)  (45) 


Consider  N2  is  inert  (N2  =  0), 
d2(lnx2)/dx2  =  Uij’/cD12  effneF.  Therefore, 


*2(*)  =  *2,bexp 


V\Ot 


c^12,effneF  L 


rj(x)-rj(0)-  —  x 

Ojnn 


that  is, 

(39) 


D 1X1  +  D2x2  =  D ix1>b  +  D2x2  b,  D3x3  +  D4X4  =  D3x3  b  +  D4x4  b 

(46) 

where  a  =  D1xlb  +D3x3  b  -  crtr](0)lcneF.  Assume  the  water-gas  shift 
reaction  is  in  thermodynamic  equilibrium  (rWG s  =  0),  i.e. 


4.4.  O2-N2-FI2O  system 


*1*4  =  Kshiftx2x3 


(47) 


It  is  a  typical  case  in  cathode  catalyst  layer  of  PEMFC  (1: 
02,  2:  N2,  3:  H20).  Assume  the  net  water  and  gas-species  flux 
via  the  membrane  is  zero  for  simplicity,  that  is,  N2  is  inert 
(N2  =  0)  and  02-H20  flux  are  related  to  as  N!/Ui=N3/u3.  So, 
d2(lnx2)/dx2=(u1/D12,eff+  ^3 1 D23teff)j I cneF,  which  leads  to 

x2(x) 


where  KShift  is  the  equilibrium  constant  of  water  shift  reaction, 
we  obtain  the  concentration  distribution  of  species  from  Eqs. 
(45)-(47), 


Xi(x)  = 


i>(x)  -  ^(x)  -  4k(k  -  1  )(D!Xljb  +  D2x2  b)/(x) 
2D,(k-  1) 


where  k  =  /CshiftDi  D4/D2D3  and 


(48) 


,  a t  v\  v3 

:X2’b6XPl  rieF  (  + 


d(x)  —  d(0)  —  — — x 

^ion  J 


(40) 


Xi(x) 


(x,. 


+  +^r T-^*2 


v\  +  v3  m  -  1 )  \  x2  b  J  v\  +  v3  m-  1 


(41) 


where 


b{x)  =  (/<-!  )/(x)  +  k(DiX1?b  +  D2x2  b)  +  D3x3?b  +  D4x4?b 


(49) 


The  concentrations  of  other  species  can  be  obtained  by  substi¬ 
tuting  x\  into  Eqs.  (45)  and  (46).  Like  the  other  cases  mentioned 
above,  the  profile  of  gas  concentrations  can  be  used  to  calculate  the 
concentration  loss  and  limiting  current  density  based  on  Fields  law. 

More  importantly,  an  analytical  expression  of  the  proportion  of 
H2  electrochemical  current  density  is  available  from  the  explicit 
profiles  of  H2  concentration  and  overpotential,  that  is, 


<JL>{X)  = 


(cDi  neF/<7t){d2X]  /  dx2) 
(d2r]/dx2) 


(50) 


(ui  +  v3 )F)23  effP]2 

D13,eff{vlD23,eff  +  ^3^12,eff)  ’ 


V\  m 

a  =  — - - 

V\  +  v3 


D13,eff\ 

^12,eff  J 
(42) 


4.5.  H2-H2O-CO-CO2-N2  system 


Instead  of  the  Stefan-Maxwell  equation,  we  use  Fields  law  for 
this  case  to  describe  the  multi-component  (1 :  H2,  2:  FI20, 3:  CO,  4: 
C02, 5:  N2)  mass  transfer  using  a  simple  analytical  solution,  that  is, 

,  5 

Ni  =  — cDj-^+Xj^JV,-  (i  =  l,  2,3,4)  (43) 

j=  1 


For  analysis  of  H2  and  CO  electrochemical  co-oxidation,  gener¬ 
ally  different  exchange  current  densities  are  used  for  FI2  and  CO, 
but  the  uncertainties  of  these  kinetic  data  are  currently  under  wide 
discussion.  On  the  other  hand,  because  the  water  gas  shift  reaction 
(WGSR)  is  fast  and  usually  assumed  to  be  in  equilibrium,  it  is  not 
possible  to  separate  the  different  exchange  current  densities  from 
an  analysis  of  gas  concentrations  only.  Still,  the  present  framework, 
resulting  in  Eq.  (50),  provides  some  information  for  kinetic  analy¬ 
sis.  Here,  the  unique  exchange  current  density  can  be  explained  as 
a  global  or  overall  one  corresponding  to  the  overall  overpotential 
driving  both  H2  and  CO  oxidation.  Fig.  5  shows  the  proportion  of  H2 
electrochemical  current  density  in  the  total  current  density,  00,  in 
the  active  zone  under  the  conditions  of  different  fuel  compositions. 
The  active  zone  is  also  called  boundary  layer  of  electrochemical 
reaction,  where  99%  variation  of  the  ionic  current  density  occurs 
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Fig.  5.  Ratio  of  H2  over  the  total  electrochemical  current  density  in  the  active  zone 
for  different  fuel  composition  of  case  1:  31.77%  H2-18.23%  H20-28.23%  CO-1 6.77% 
C02-5%  N2>  J=2  x  104  Am-2,  case  2:  22.24%  H2-7.76%  H20-47.76%  CO-17.24% 
C02-5%  N2> /=  1  x  104  A rrr2,  case 3:  8.98%  H2-1.02%  H20-76.02%  CO-8.98%  C02-5% 
N2, 7  =  0.5  x  104  Am"2,  and  the  other  parameters  are  the  same  as  those  in  Fig.  2. 


[17].  It  is  obvious  that  the  fuel  composition  has  a  great  influence 
on  the  electrochemical  kinetics  of  H2  and  CO.  As  CO  concentration 
increases  or  H2  concentration  decreases,  oo  decreases.  However, 
when  the  concentrations  of  H2  and  CO  are  comparable,  most  of 
electrochemical  current  density  (over  80%  in  case  1)  results  from 
oxidation  of  H2. 


The  heat  source  is  given  by 


Qk  —  Sh.khk(Ts  Tk)  +  Smk 


^  ^ ref^i,  a^ref  +  ^i,WGS^i,a^WGs) 


'y  "max(  ^rib^i,E/C5  ^)^p,iiTs  T) 

i  J  k 


(54) 


where  Ta,  Tc,  Ts  are  temperature  of  the  anode  gas,  cathode  gas  and 
solid  matrix,  u  is  the  gas  velocity,  h  is  the  convective  heat  transfer 
coefficient,  x,-,  Cp  i  and  are  the  molar  fraction,  molar  specific  heat 
and  enthalpy  of  species  i,  Cp  =  ^  iCPtiXi  is  the  gas  molar  specific  heat. 
For  rectangular  flow  channels  in  a  typical  planar  SOFC,  the  specific 
mass  transfer  area  (Sm)  and  heat  transfer  area  (Sh)  per  unit  volume 
of  flow  channel  and  rib  coefficient  (§rib)  are 


Sm,k  — 


Jch,k 


^h,k  — 


2(DCh,k  +  Wch) 
Dch,kWch 


£rib  —  1  + 


Wr 


rib 


w( 


ch 


(55) 


where  Dch,  Wch  are  the  depth  and  width  of  flow  channel  and  Wrib 
is  the  rib  width. 

For  the  solid-phase  energy  balance  it  is  assumed  that  tem¬ 
perature  varies  only  in  z  direction,  while  in  x  direction  all  solid 
components  (MEA  and  interconnectors)  have  the  same  tempera¬ 
ture.  Radiant  heat  transfer  is  neglected.  The  energy  balance  is  given 
by 


rO.  t 

*s  ^2  =  ^  y  Kk^rib^m.k  E/C^k,  i  ( ^k  ?  Ts . 

k=a,c  i 


5.  Along-the-channel  cell-level  model 

The  approximate  analytical  solution  of  electrode-level  model 
derived  above  is  further  integrated  into  a  ID  thermal  cell-level 
model  where  transport  of  mass  and  heat  along  a  single  channel  is 
considered.  This  yields  an  overall  1 D  + 1 D  model,  where  the  dimen¬ 
sion  x  through  the  MEA  is  solved  analytically  and  the  dimension  z 
along  the  channel  length  is  solved  numerically.  Only  steady  state  is 
considered  here. 

5.1.  Governing  equations  and  background 

Along  the  coordinate  of  channel  length,  z  e  [0,  L],  the  gas-phase 
energy  balance,  species  and  overall  mass  balance  are  given  by  [4] 

Ck«kCp,k^=QK  ( k  =  a,c )  (51) 

=  5m,k(— £rib^i,E/C  3"  rcf^ref  3"  WCS^WGs)^  (52) 

=  ^k5m,k^~7~£ribMj,E/C  3“  ref^rcf  3“  WGS^WCS  )^  9-  ~q  ~ 

(53) 


Here,  Eq.  (53)  is  a  combination  of  the  total  mass  balance  and  the 
ideal  gas  law,  p  =  cRT ,  so  dp/dt  =  R(Tdc/dt  +  cdT/dt ).  Substituting  Eqs. 
(51 )  and  (52)  to  this  expression,  one  will  find  that  there  is  no  tran¬ 
sient  term  (even  for  transient  simulation)  when  not  considering  the 
pressure  dynamics  (3p/3t  =  0,  dp/dz  is  not  forcedly  excluded).  The 
advantage  of  this  formulation  is  to  provide  an  equation  with  gas 
velocity  included  as  dependent  variable,  as  we  do  not  put  forward 
the  complex  momentum  equation  (like  in  the  full  Navier-Stokes 
equations).  This  treatment  has  been  shown  before  to  be  numeri¬ 
cally  much  more  stable  than  the  direct  use  of  the  ideal  gas  law  [4], 
and  can  be  also  found  in  Heidebrecht’s  work  [25]. 


y  ^  Kk^h.k^kC^ik  rsJ  +  ycfribWWcdl  (56) 

k=a,c 

where /cs  is  heat  conductivity,  yk  =  Dch)kWch/(AC0  n,a  + A:on,c  +  ImeaW ch) 
the  ratio  of  solid-phase  heat  conduction  area  to  gas-phase  mass 
or  heat  transfer  area,  where  /MEA  is  the  MEA  thickness,  Acon  the 
section  area  of  interconnector. 

The  cell  voltage,  Vcen,  which  is  assumed  constant  along  the  chan¬ 
nel  (ideal  conductivity  of  the  interconnectors),  can  be  obtained  by 

VceH  =  Voc (Z)  -  l/ta(z)  -  l/tc(z)  -  ^  (57) 

°ion 

where  /(z)  is  the  local  cell-area-specific  current  density,  le  the  thick¬ 
ness  of  the  electrolyte  layer,  the  local  open  circuit  voltage,  V0c  is 
calculated  by  the  Nernst  equation 


£o,ref  + 


^5ref 

neF 


(Ts 


Tref )  + 


RT  ^Xm,b(PcXQ2,b/Po)°'5 
neF  *H20,b 


(58) 


where  E0,ref>  ASref  are  the  reference  electromotive  force  and 
entropy  change  of  the  overall  electrochemical  reaction  (here: 
H2  +  0.502  ->  H20),  and  the  total  overpotential  everywhere  along 
the  channel,  rjt(z)  can  be  obtained  from  the  approximation  solution 
p(z,x)  and  Xj(z,x)  of  the  local  electrode-level  model  via 


rh(z)  =  — E— [ffion ri(z,  0)  +  <Te\tj(z,  l)  +  I(z)l]  +  ^ln^^-(59) 

The  overall  equation  system  represents  the  implicit  relation¬ 
ship  between  cell  voltage  Vcen  and  average  cell  current  density  Jcell. 
Either  can  be  independently  specified,  and  the  other  value  is  sim¬ 
ulated.  The  fuel  cell  is  generally  treated  as  a  set  of  parallel  discrete 
elements  with  uniform  cell  voltage  and  non-uniform  local  current 
densities.  Thus,  as  shown  in  Fig.  6,  a  high-level  iteration  loop  is 
required  for  the  cell  voltage  or  average  current  density.  This  itera¬ 
tion  loop  greatly  increases  the  computational  cost  and  also  the  risk 
of  numerical  divergence,  although  this  problem  can  be  alleviated 
by  equation-oriented  solvers  like  gPROMS  [4]. 
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Fig.  6.  Flow  chart  of  voltage-current  iteration  loop  in  distributed  fuel  cell  modeling. 


5.2.  Discrete  form  and  model  simplification 


xH2,j  =  XH2  J-l  +  ,  *02, j  = 


(x02  J— 1  +  ^02 {/—l  ) 
(1  +  /<02^j-l) 


Assuming  that  the  total  gas  pressure,  gas  molar  specific  heat,  and 
electrochemical  reaction  heat  (AH)  are  constant,  we  further  sim¬ 
plify  the  computational  framework  of  distributed  modeling.  For  the 
simplest  case  in  SOFCs,  that  is,  FI2/FI2O  at  the  anode  and  02/N2  at 
the  cathode,  firef  =  KWgs  =  0,  %/c  =  -ViI{z)lneF,  caua  =  caiinuaiin,  and 
the  variation  of  total  gas  molar  flux  in  the  cathode  flow  channel  is 
considered  negligible  (ccuc  =  cc  inuc  in)  due  to  a  generally  large  stoi¬ 
chiometry  of  02  and  dilute  effect  of  N2 .  Based  on  the  control-volume 
discretization  approach  as  shown  in  Fig.  7,  we  can  obtain  the  fol¬ 
lowing  discrete  forms  according  to  the  power  law  or  the  upwind 
scheme  [26]. 

For  co-flow, 

Ta,j  ~  tzTaj_i  +  (1  -  o)Tsj,  Tcj  =  bTcj_i  +  (1  -  b)TSj 
0  =  2,  ...,n  +  l)  (60) 


and  for  counter-flow, 

■Ta  ,j  —  a^a,j-i  +  (1  —  a)TsJ ,  Tcj  =  bTcj_ |_|  +  (1  —  b)Tsj 
0  =  2,  n  +  1) 


*H2  ,j  =  *H2 J-l  +  ^H2fj-1 , 


_  (*02  ,j+l  +  ^02^-1  ) 

02J'  “  (1+V/-1) 


where 
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(1  +  (^h.c^^c/Ccjn^c.in^p.c)) 
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Fig.  7.  Discrete  scheme  of  control  volumes  used  in  the  present  work. 
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Fig.  8.  Comparison  between  pressure-dependent  V-I  performances  of  the  approxi¬ 
mation  solution  (this  work)  and  numerical  elementary  kinetic  simulation  [27].  The 
parameters  are:  k  =  540  pm,  lc  =  55  pan,  le  =  10  pm,  Wch  =  2  mm,  Dcha  =Dchc  =  2  mm, 
Wrib  =  2  mm,  rPia  =  0.5fjim,  rPx  =  0.25  pm,  ea  =  0.32,  ra  =  2.9,  ec  =  0.4,  rc  =  1.4, 
orion  =  5.15  x  107/T exp( -84, 000 /RT) S m-1 ,  <xeU  =  9.5  x  107/Texp(-1150/T)Sm-1, 
orei,c  =  4.2x  107/Texp(-1200/T)Sm-1,  £act,a  =  130kjmol1,  Fact>c  =  125kJmol-1, 
o(a  =  ac- 0.5,  io,H2  =  1  x  108  Am ~3,  Vh2=0.5,  Vh2o  =  0»  T=Tain  =  TCin  =  rref=1073.15K, 
feq.fuei  =  2  x  1 04  A  m-2 ,  Jeq  air  =  4  x  1 04  A  m-2 ,  and  the  cathode  overpotential  is 
calculated  by  interfacial  BV  model  as  /o,o2  =3.6  x  104  Am-2,  v02  =  0.466. 


Fig.  9.  Comparison  between  V-I  performances  of  the  approximate  simulation  (this 
work)  and  a  full  numerical  model  [4]  under  different  fuel  compositions  for  case  1: 
97%  H2-3%  H20,  case  2: 48.5%  H2-1.5%  H20-50%  N2,  case  3: 19.4%  H2-0.6%  H2O-80% 
N2. 
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Fig.  10.  Comparison  between  1D  +  1D  simulation  and  ID  simulations  (through  electrode  thickness).  The  ID  model  uses  as  boundary  conditions  for  bulk  concentration  at 
channel/electrode  interface  (a)  ID  +  ID  inlet  species  concentrations,  (b)  ID  +  ID  outlet  species  concentration,  (c)  arithmetic  mean  and  (d)  logarithmic  mean  between  1D  +  ID 
inlet  and  outlet  species  concentration.  Parameters  are  identical  to  those  in  Fig.  8  and  the  fuel  compositions  are  identical  to  the  three  cases  in  Fig.  9. 
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Fig.  11.  Along-the-channel  profiles  of  local  current  density,  solid  and  gas  temperatures  and  overpotentials  of  both  exact  and  approximate  calculations  in  the  case  of  (a) 
concurrent  flow,  (b)  countercurrent  flow.  Fuel  composition  is  97%  H2-3%  H20  and  air  is  used  as  oxidant.  For  concurrent  flow,  7ain  =  Tcin  =  975.15  K,  and  for  countercurrent 
flow  ra  in  =  Tc jn  =  958.15  K,  other  basic  parameters  are  identical  to  those  in  Fig.  8. 


Substituting  Eq.  (60)  or  (62)  into  Eq.  (56),  we  can  get  a  lin¬ 
ear  relationship  between  the  local  current  densities  and  solid 
temperatures,  [I]  =  [A][TS]  +  [B],  which  is  explained  in  Appendix  B. 
In  this  case,  the  calculating  mode  with  given  Vceli  (cf.  Fig.  6)  is 


adopted.  Substituting  the  gas  concentrations  (Eq.  (61)  or  (63))  as 
functions  of  Ts(z)  or /(z)  into  Eqs.  (57)-(59),  a  fast  and  stable  along- 
the-channel  computation  can  be  obtained  without  V-I  iteration 
loop. 
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6.  Results  and  discussion 

6.1.  Isothermal  ID  +  ID  simulations 

We  first  compare  the  present  model,  which  is  based  on  a 
number  of  approximation  assumptions,  to  two  different  more  com¬ 
plex  literature  models.  Fig.  8  presents  the  pressure-dependent  V-l 
performance  of  an  anode-supported  SOFC.  Our  approximate  sim¬ 
ulation  is  compared  to  a  model  by  Henke  et  al.  [27],  which  is 
implemented  in  an  in-house  C-code  software  package  DENIS  for 
detailed  electrochemistry  and  numerical  impedance  simulation. 
All  the  cases  are  isothermal  and  under  countercurrent  flow.  The 
inlet  fuel  and  oxidant  flow  rates  are  set  to  a  current-density  equiv¬ 
alent  (Jeq)  of  2  and  4  A  cm-2,  respectively.  The  present  model  was 
adjusted  to  the  numerical  model  by  fitting  the  volumetric  exchange 
current  density  in  the  anode,  anode  average  particle  size,  and  H2 
reaction  order.  The  overestimation  of  the  present  model  in  the  acti¬ 
vation  zone  is  probably  due  to  the  different  electrochemical  kinetics 
frameworks  (one  is  based  on  overall  Butler-Volmer  equation  and 
the  other  is  based  on  elementary  kinetics).  On  a  2.26  GHz,  1.92  GB 
PC,  the  computational  time  is  0.87  s  per  point  along  the  IV-curve 
for  the  present  Matlab-code  approximate  model  (20  along-the- 
channel  control  volumes)  and  average  430  s  per  point  (10  control 
volumes  along  the  channel  length  and  22  control  volumes  along  the 
anode  thickness)  for  the  detailed  Henke  et  al.  model.  These  numbers 
demonstrate  the  high  performance  of  the  present  algorithm. 

Fig.  9  further  compares  the  V-I  performance  between  the 
approximate  simulation  and  numeric  calculation  by  a  model  of  Bao 
et  al.  [4].  The  latter  is  based  on  a  gPROMS  implementation  of  a  multi¬ 
level  SOFC  modeling  platform  [4],  which  uses  the  identical  physics 
as  described  in  the  present  paper,  however  without  approximating 
assumptions.  The  comparison  is  carried  out  for  three  different  fuel 
compositions.  All  the  cases  are  isothermal  and  under  countercur¬ 
rent  flow,  and  the  cell  configuration  and  simulation  parameters  are 
identical  to  those  used  in  Fig.  8.  The  two  simulations  are  almost 
identical  as  predicted  both  under  open  circuit  conditions  as  well 
as  under  diffusion-limited  conditions.  The  Bao  et  al.  model.  [4] 
(20  along-the-channel  control  volumes  and  15  along-the-anode 
control  volumes  based  on  orthogonal  collocation  finite  elements 
discretization  method)  takes  average  15  s  for  each  point  along  I-V 
curve.  Therefore,  the  approximate  approach  is  dozens  of  times 
faster  than  the  equation-based  solver  in  gPROMS  and  furthermore 
found  to  be  more  robust  in  the  activation  and  concentration  loss 
zone. 

6.2.  1 D  simulations  using  effective  bulk  concentrations 

As  exemplary  application  of  the  present  model,  we  compare  ID 
with  ID  +  ID  isothermal  simulations.  For  a  potential  further  model 
reduction,  it  is  interesting  to  investigate  under  which  boundary 
conditions  a  simple  ID  model  can  best  represent  a  more  com¬ 
plex  1D  +  1D  simulation.  In  particular,  we  compare  four  different 
options  for  the  species  bulk  concentration  at  the  channel/electrode 
interface,  (i)  using  the  inlet  concentration,  (ii)  using  the  outlet  con¬ 
centration,  (iii)  using  the  arithmetic  mean  between  inlet  and  outlet 
concentrations,  and  (iv)  using  the  logarithmic  mean  between  the 
inlet  and  outlet  concentrations.  Fig.  10  compares  the  1D  +  ID  sim¬ 
ulation  with  these  four  kinds  of  ID  simulations  under  different 
fuel  compositions.  As  Fig.  10(a)  shows,  taking  the  inlet  concentra¬ 
tion  leads  to  a  considerable  overestimation  of  the  limiting  current 
density.  In  this  case,  consumption  of  fuel  and  oxidant  cannot  be 
reflected  reasonably,  which  results  in  the  species  concentration 
at  the  electrode/electrolyte  interface  even  higher  than  the  outlet 
concentration.  According  to  the  upwind  scheme  in  1D  +  1D  mod¬ 
eling,  the  outlet  species  concentration  is  taken  as  the  bulk  one 
in  each  control  volume,  i.e.  each  element  cell  is  considered  as  a 


continuous  stirred-tank  reactor  (CSTR).  However,  Fig.  10(b)  shows 
that,  using  the  outlet  concentration  in  ID  simulation  results  in  an 
obvious  underprediction  of  cell  performance.  The  method  of  arith¬ 
metic  mean  provides  a  balance  between  the  first  two  options,  but 
still  leads  to  an  overestimation  of  the  limiting  current  density  as 
shown  in  Fig.  10(c).  As  shown  in  Fig.  10(d),  the  logarithmic  mean 
seems  to  be  the  best  one  with  regard  to  the  limiting  current  den¬ 
sity.  The  physical  meaning  of  logarithmic  mean  concentration  can 
be  compared  to  the  logarithmic  mean  temperature  difference  in 
heat  exchangers.  Bao  et  al.  have  previously  used  the  concept  of  log¬ 
arithmic  mean  bulk  concentrations  in  modeling  of  PEM  fuel  cells 
and  systems  [28,29]. 


6.3.  Full  thermal  1D  +  1D  simulations 

Finally,  we  apply  the  full  thermal  1D  +  ID  model.  Fig.  11  shows 
along-the-channel  profiles  of  the  local  current  density,  solid  and 
gas  temperatures  and  overpotentials.  The  results  of  the  approxi¬ 
mate  model  (this  work)  are  compared  to  exact  solutions  (Bao  et  al. 
model  [4])  under  non-isothermal  condition.  For  these  simulations, 
the  inlet  fuel  and  air  temperatures  were  adjusted  to  values  of  702 
and  685  °C,  which  makes  the  average  cell  temperature  approxi¬ 
mately  800  °C.  Here,  the  exact  solution  in  gPROMS  [4]  considers  the 
radiant  heat  transfer  between  the  MEA  and  interconnectors  with 
detailed  view  factor  models,  while  it  is  excluded  in  the  approximate 
calculation  due  to  the  assumption  of  identical  local  temperature  in 
the  both  solid  phases.  As  shown  in  Fig.  11(a)  and  (b),  the  approxi¬ 
mate  model  as  derived  in  the  present  work  agrees  well  with  the 
exact  solution  under  both  co-flow  and  counter-flow  conditions. 
Moreover,  the  approximation  method  is  found  to  be  rather  sta¬ 
ble  even  with  a  coarse  mesh  grid.  The  present  model  with  20 
or  50  along-the-channel  control  volumes  takes  only  1.3  and  6.7  s 
respectively  to  complete  a  full  thermal  1D  +  ID  numeration,  while 
it  correspondingly  takes  average  1 00  or  340  s  for  the  gPROMS-code 
exact  computation. 


7.  Conclusions 

As  balance  between  full  mechanistic  and  semi-empirical  fuel 
cell  models,  approximate  analytical  solutions  provide  a  computa¬ 
tionally  efficient  framework  with  a  clear  physical  meaning.  This 
is  an  important  aspect  for  multi-scale  models  and/or  real-time 
simulations  of  fuel  cells,  where  a  compromise  needs  to  be  found 
between  computational  speed  and  model  precision.  In  the  present 
work,  analytical  solutions  for  ID  through-the-electrode  charge 
and  mass  transport  were  developed  and  integrated  into  a  ID 
along-the-channel  thermal  numerical  simulation.  The  approxima¬ 
tion  approaches  shown  in  this  paper  keep  the  intrinsic  nonlinearity 
of  system. 

For  charge  transport  and  reaction  inside  a  thin  porous  electrode, 
overpotential  shows  a  smooth  profile  and  perturbation  method 
generally  fails  to  find  a  small  enough  variable.  Therefore,  a  mod¬ 
ified  Adomian  decomposition  method  (ADM)  was  applied  and 
validated.  Furthermore,  an  explicit  hybrid  algorithm  of  power-law 
approach  and  linear  approximation  based  on  WKB  perturbation 
method  was  introduced  to  overcome  some  disadvantages  in  the 
ADM  solution  (too  many  terms,  lacking  of  clear  physical  mean¬ 
ing,  implicit  iteration).  Although  starting  from  a  thick  electrode, 
this  algorithm  is  validated  to  accurately  reflect  the  boundary-layer 
effects  at  both  electrode  interfaces  in  a  wide  range  of  applica¬ 
tions.  In  addition,  a  novel  expression  was  developed  to  interconvert 
the  volume-specific  and  area-specific  exchange  current  densities, 
which  provides  a  better  theoretical  explanation  over  the  common 
linear  relationship  via  a  simple  geometric  factor. 
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Exact  solutions  of  mass  transfer  of  different  fuel-cell-relevant 
gas  mixtures  were  then  obtained  by  combining  Stephan-Maxwell 
diffusion  equations  with  the  overpotential  profiles  resulting  from 
the  charge  transport  model.  In  the  case  of  the  H2-H2O-CO-CO2-N2 
system  (SOFC  anode),  by  considering  H2/CO  co-oxidation,  we 
derived  an  analytical  expression  of  the  ratio  of  H2  and  CO  elec¬ 
trochemical  current  density.  In  this  regard,  the  thermodynamic 
framework  can  provide  a  first  information  for  kinetics  analysis. 

We  further  integrated  the  electrode-level  model  into  an  along- 
the-channel  simulation.  This  results  in  an  overall  1D  +  1D  model, 
where  one  dimension  is  solved  analytically  and  one  dimension  is 
solved  numerically.  For  the  reactant  system  of  H2-H20-N2  and 
02-N2,  a  linear  relationship  between  the  along-the-channel  pro¬ 
files  of  local  current  density  and  solid  temperature  was  further 
developed  in  both  cases  of  concurrent  and  countercurrent  flow. 
The  1D  +  1D  model  was  validated  by  comparison  to  two  different 
numerical  models  developed  previously  by  the  authors.  In  addition, 
when  comparing  1 D  + 1 D  against  1 D  models,  we  found  that  the  log¬ 
arithmic  mean  of  the  inlet  and  outlet  species  concentrations  can  be 
used  as  boundary  condition  at  the  channel/electrode  interface  in 
ID  simulations. 

Based  on  the  approximations  applied  to  both  dimensions,  the 
computational  efficiency  of  the  1D  +  1D  model  presented  in  this 
paper  approaches  that  of  a  zero-dimensional  model.  Generally,  it 
only  takes  several  seconds  to  complete  a  distributed  cell-level  sim¬ 
ulation,  which  is  valuable  for  real-time  simulations  and  multi-scale 
fuel  cell  models. 
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Appendix  A.  WKB  perturbation  method 

Expand  the  exponent  function  as  Taylor  series 

00  n 

/  =  fco[exp(aay)  -  exp(-acy)]  =  kp'Y'f3 

n=0 


00 

=  X2^a„y"  (Al) 

n= 1 

where 

x-v^TT^),  m 

When  k0  is  a  large  number,  1  /X  can  be  used  as  a  perturbation 
variable.  Assume 

y  =  exp(XG)  and  G  =  G0  +  A_1Gi  (A3) 

When  only  zero-order  term  of  exponent  function  is  considered 

00  00  00 

(A4) 

n= 1  n= 1  n= 1 

Substitute  Eqs.  (A3)  and  (A4)  into  Eq.  (Al ),  we  obtain 

00 

G'l  +  X-\2G'0G\  +  C")  =  (A5) 

n= 1 


5>y"  =  ^aneAnG  «  eXG^an 


So, 


G'o  - 

n= 1 

2G()G/1  +  G"  =  0 


J0  ' 


(A6) 


Gi  —  2 111  ^0 


Substitute  G0  and  Gi  into  Eq.  (A3),  we  can  obtain  the  general 
solution 


y  =  ci  exp 


I  eaa  —  e_“c 


Oi  a  +  OLq 


x  +  c2  exp  -X 


g<Xa  — 

Ctfa  +  Oi  c 


x  (A7) 


where  G\  has  been  included  in  the  constant  C\  and  c2. 


Appendix  B.  Linear  relationship  between  current  density 
and  solid  temperature 

From  Eqs.  (60)  and  (62),  there  is 


m=j- 1 

Taj  -  TSJ  =  o'-'Ta,,  +  (1  -  q)y>-rs,m  -  tTTsj 

m=2 
m=j- 1 

Tcj  -  Tsj  =  W-'Ic,,  +  (1  -  -  bTSJ  (coflow) 


(Bl) 


Tcj  -  Tsj  =  b"-J+2rc,n+2  +  (1  -  iTs'm  ~  bTsJ  (counterfIow) 


Substitute  it  into  the  discrete  form  of  Eq.  (56),  we  obtain 
the  linear  relationship  between  the  local  current  den¬ 
sity  vector  [I]  =  [Ji,. .  ./n]r  and  the  solid  temperature  vector 
[Ts]  =  [TS2,. .  •Tsn+i  i.e.  [I]  =  [A][TS]  +  [B],  the  elements  of  matrix 
[A]  and  [B]  are 

For  concurrent  flow 


A\,\  = 
M,2  = 
= 

A',i-1 

\i  = 

A,i+ 1 
Anj  = 
An,n-\ 
An,n  = 


-{YaSh,ahaa  +  ycShiChcb  +  (/cs/Az2)) 

^  ,c 
KAz 2  ’ 

[yaSh.aha(l  -  a)a'-J+2  +  yc5h,cftc(l  -  b)bl~^2] 

_  [XaSh,a/ia(l  -  a)a  +  ycSh,cA(  1  -  b)b  +  (/Cs/Az2)] 
K 

-(KaSh,a/iaa  +  YcShjChcb  +  (2/cs/Az2)) 


(B2) 


KAz 2 


(i  =  2 . .  .n  -  1) 


[yaSh,aha(l  -  a)an~j+2  +  ycSh>c/ic(l  -  b)bn~>+2] 


I< 


(j=  1...H-2) 


lYaShMl  -  a)a  +  XcSh,chc(l  -  b)b  +  (/cs/Az2)] 
I< 

-(XaSh,a/Taa  +  ycShxhcb  +  (/cs/Az2)) 

I< 


D  (Ka^h.a^a^T’a.in  +  Yc^h,c^cb^Tc[n) 

B,= - K -  0  =  1 ' 


■  •  n) 


(B3) 
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For  countercurrent  flow 

— (Xa^a^a#  +  Yc$h,cbcb  +  {ks/ Az2)) 


A,1  = 


[ycSh,A(l  —  +  (/rs/Az2)] 


/C. 


...-VAAa-tife-'  tf,3  n) 

0  —  1...,  — 2) 

*  [yaSh,afia(l -a)a  +  (/cs/Az2)]  0  „  ,, 

A,i-i  = - ^ -  (i  =  z...n  i j 


/< 

— (/a^h.a^aQ  +  ycShChcb  +  (2/cs/Az  )) 
[XcSh.cWl  —  b)b  +  (/Cs/Az2)] 

tf,i+2...n) 


Am  - 


A,t+i  = 


(B4) 


Vn-1  = 


Ai,n  = 


/C 

[ya5h,a^a(l  —  +  XcSh.cfrcO  —  b)b  +  (/Cs/Az2)] 

k  : 

— (Ka-^h.a^aQ  +  Kc-^h.c^c^  +  (^s/Az  )) 

K 


D  (ya5h,a^aaI7a,in  +  Yc^h,cbcbn  I+^c,in)  ^  ^ 

^ -  U  =  1  •  •  • 


(B5) 


where  the  adiabatic  conditions  of  solid  phase,  3Ts/dz|z_0  = 
3Ts/9z|2_l  =  0  have  been  included  and  the  denominator  JC  is 

K  =  ■WcSm.c  (VCell  +  j  (B6) 
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